Design and analysis of a flexible Ruddlesden–Popper 2D perovskite metastructure based on symmetry-protected THz-bound states in the continuum

A Ruddlesden–Popper 2D perovskite PEA2PbX4 (X = I, Br, and Cl) is proposed for metasurface applications. Density functional theory is used to analyze the optical, electrical, mechanical properties, moisture and thermodynamic stability of PEA2PbX4. The refractive index of PEA2PbX4 varies with the halides, resulting in 2.131, 1.901, and 1.842 for X = I, Br, and Cl, respectively. Mechanical properties with Voigt-Reuss-Hill approximations indicate that all three materials are flexible and ductile. Based on the calculations of formation energy and adsorption of water molecules, PEA2PbI4 has superior thermodynamic and moisture stability. We present a novel metasurface based on 2D-PEA2PbI4 and analyze symmetry protected-bound states in the continuum (sp-BIC) excitation. The proposed structure can excite multiple Fano quasi-BICs (q-BICs) with exceptionally high Q-factors. We verify the group theoretical analysis and explore the near-field distribution and far-field scattering of q-BICs. The findings indicate that x-polarized incident waves can excite magnetic toroidal dipole-electromagnetic-induced transparency-BIC and magnetic quadrupole-BIC, while y-polarized incident waves can excite electric toroidal dipole-BIC and electric quadrupole-BIC. The influence of meta-atom and substrate losses, array size limitations, and fabrication tolerances are also discussed. The proposed structure can be employed for applications in the THz region, such as polarization-dependent filters, bidirectional optical switches, and wearable photonic devices.

The 2D perovskites exhibit unique and remarkable structural and physical properties, such as low cost, ease of fabrication, direct and tunable bandgap, a soft and dynamic structure, and a relatively high nonlinear refractive index.Unlike 3D perovskites, the flexibility of the organic layer in 2D perovskites further contributes to their softness, along with the Van der Waals interface.In addition, 2D perovskite materials exhibit higher environmental and chemical stability compared to their 3D counterparts 26 .In the 2D perovskite material PEA 2 PbX 4 (X = I, Br, and Cl), Phenethylammonium (PEA) with the chemical formula C 6 H 5 (CH 2 ) 2 NH 3 is used as the large cation.These aromatic cations exhibit hydrophobic characteristics typical of bulky aromatic cations.Consequently, they not only improve the material's environmental stability but also can influence its mechanical properties.In layered 2D halide perovskites based on PEA, moisture stability is generated by the hydrophobic nature of aromatic organic ammonium spacer cations.In general, perovskite materials, especially 2D perovskite (i.e., PEA 2 PbX 4 ), have more flexibility than silicon because they have smaller elastic moduli 27,28 .PEA 2 PbX 4 has recently been used in many optoelectronic fields due to its outstanding properties, such as moisture stability, photostability, and π-π interaction 29 , and we want to introduce and study it in the THz region for the first time to the best of our knowledge.
In the present work, we employ first-principles density functional theory (DFT) analysis to examine the mechanical, optical, electrical, and stability characteristics of PEA 2 PbX 4 (X = I, Br, and Cl) to investigate their potential use in THz applications.The results indicate that PEA 2 PbI 4 is significantly more thermodynamically stable than PEA 2 PbBr 4 and PEA 2 PbCl 4 , and is some orders of magnitude more stable than 3D perovskites.The DFT calculation is used to assess the shear modulus, Poisson's ratio, Young's modulus, and bulk modulus of these compounds.The elastic moduli calculated by the Voigt-Reuss-Hill approximations indicate that these compounds have both ductility and mechanical stability.Moreover, for the first time, we demonstrate the moisture stability of PEA 2 PbX 4 by calculating the adsorption energy of water molecules.We find that these materials have a high dielectric constant, near-zero loss in the THz region, thermodynamic stability, moisture resistance, and flexibility.
Following that, we introduce a novel kind of metastructure that uses periodic circular slot rings in a Ruddlesden-Popper quasi-2D perovskite layer (PEA 2 PbI 4 ), which has not been reported yet to the best of our knowledge.By introducing an offset distance of the inner ring from the center, the symmetry of the structure is disrupted, which opens a zero-order channel, enabling the conversion of the dark modes to q-BIC states with a finite linewidth.After conducting an extensive mathematical examination utilizing group theory and finite element eigenfrequency calculations, we prove that the proposed metastructure can excite multiband high-quality q-BICs, such as the magnetic toroidal dipole (MTD) and electric toroidal dipole (ETD), which exhibit specific symmetry properties.By calculating the far-field radiations from the multipole decomposition method and nearfield analysis with finite element frequency domain calculations, a thorough study is conducted on the q-BIC characteristics.The analysis of symmetry is beneficial in determining which BICs can be excited based on the polarization of the incident wave.Moreover, when the broad mode is coupled to the q-BIC, the Fano peak of low bandwidth transparency is visible.The Q value of electromagnetic-induced transparency-BIC resonance can be regulated by adjusting the offset distance.This work presents a valuable reference for developing applications in the THz region such as polarization-dependent switches, multi-channel wearable biochemical sensing, notch polarization-dependent filters, and non-linear optics.

Materials and methods
The schematic crystal structure of PEA 2 PbX 4 (X = I, Br, and Cl) is depicted in Fig. 1a.Table S1 presents crystal data for room temperature structures.

Optical properties
The first-principle analysis based on DFT is used to compute the complex refractive index (RI) of PEA 2 PbX 4 (X = I, Br, and Cl).The dielectric function, denoted as ε(ω) = ε 1 (ω) + iε 2 (ω), can be expressed as a complex quantity with a real part ε 1 (ω) and an imaginary part ε 2 (ω).To obtain optical properties, the imaginary part of the dielectric function is computed using the momentum matrix elements that correspond to the unoccupied and occupied wave functions in agreement with the selection rules, and is expressed as 30 : where ħω is the photon energy, e represents the electronic charge, ε 0 is the dielectric constant in vacuum, Ω is the volume of a unit cell, u is the unit vector along the polarization of the incident electric field, and ψ k v and ψ k c are wave functions for valence and conduction band electrons at a certain wave number, respectively.By employing the Kramers-Kronig relations, it is possible to extract the real component of the dielectric function from the imaginary component.The relation between complex dielectric function and complex RI is 30 : where k and n are the imaginary and real parts of the refractive index, respectively (see Sect.S1).

Thermodynamic stability
The formation energy (FE) of quasi-2D perovskite ( L 2 A n−1 B n X 3n+1 ) is defined as 31 : where E is the total energy of the component material.For pure 2D perovskites (n = 1), the equation is modified as: More negative formation energies directly correlate with a higher level of thermodynamic stability in the system 31 .

Mechanical properties
The mechanical properties of perovskites are critical for the analysis of ductility and flexibility in various applications.The mechanical properties are calculated using the stress-strain method to compute the elastic constant elements (C ij ).From Hooke's law, the stress (σ i ) and strain (ϵ j ) in Voigt notation are related as 32 : The shear modulus (G), bulk modulus (B), Young's modulus (E), and Poisson's ratio (τ) are calculated by DFT analysis.Employing Voigt (B V , G V , E V , ν V ) and Reuss (B R , G R , E R , ν R ) approximations, we can define the relations of elastic moduli by [33][34][35][36] : (1) Pugh's ratio refers to the ratio of the bulk modulus to the shear modulus (B/G).If the ratio of B/G is above 1.75, the material is characterized as ductile; otherwise, it is considered brittle.Another parameter for separating ductile from brittle materials is Poisson's ratio.Ductility is observed in materials with a Poisson's ratio greater than 0.26, whereas materials with a Poisson's ratio of 0.26 or lower exhibit brittleness.Generally, materials with a higher Poisson's ratio are more ductile 37 .

Moisture stability
The adsorption energy of water on perovskites is calculated by 38 : where E adsorbate presents the total energy of adsorbate (i.e., water molecules), E sub is the total energy of the isolated substrate system (i.e., perovskite surface), and E adsorbate/sub is the total energy of adsorption system (i.e., water on the perovskite surface).For the negative value of E ads , perovskite has hydrophilicity.On the other hand, for a positive E ads value, perovskite has hydrophobicity.

Structure design and symmetry analysis
The proposed metastructure comprises a PDMS substrate (n = 1.4) with a thin layer of Ruddlesden-Popper perovskite on top, which has a periodic arrangement of slot rings (Fig. 1a).The symmetric and asymmetric unit cells (UC) are illustrated in Fig. 1b and c, respectively, and are characterized by certain structural features such as the inner and outer radii of the rings (r 1 = 67.5 µm and r 2 = 75 µm), the gap between adjacent rings (g = 31.25 µm), the thickness of the perovskite layer (h = 72.5 µm), and the period of the UC (P = 181.25 µm).By adjusting the inner ring position, an asymmetric structure is achieved.The offset distance of the inner ring from the center is denoted by d.Typically, the lithography process affects mainly the roughness and uniformity of the structure while the position of the rings remains fixed.This allows precise control over the asymmetry parameter during experimentation.This approach does not change the volume of the material portion of the metasurface, resulting in a minor shift in resonance frequency that is mainly caused by variations in the coupling between neighboring rings.Overall, this method results in a relatively stable resonance frequency for q-BICs.To determine the optical characteristics of the metastructure, the finite element method is utilized.Floquet-Bloch periodic boundary conditions are used along the x-z and y-z planes.Additionally, two perfectly matched layers are positioned at a specified distance from the structure along the z-axis and are supported by scattering boundary conditions.
The unperturbed structure's UC is depicted in Fig. 1d, displaying a 2D group of geometrical symmetry denoted as C 4v in Schoenfies notation 39 .Table 1 shows the four dark modes with irreducible representations (IRREPs) A 1 , A 2 , B 1 , and B 2 , and two degenerate bright modes (E), which are supported by the C 4v symmetric group.Table 2 displays the resonant frequencies and field profiles associated with each mode.It should be noted that we use PEA 2 PbI 4 in eigenfrequency calculations; the reason for this will be explained in the next sections.Four dark modes (i.e., real eigenfrequency) have different symmetries compared to the incident electric field and cannot be excited by it, as they belong to different IRREPs 40 .Because their components are orthogonal to that of the electric field and cannot be activated by incident waves.The arrow plots presented in Table 2 provide a visual representation indicating that these dark modes lack a total electric dipole moment in the x-y plane, thereby making them incapable of being coupled with the excited plane wave.Additionally, IRREPs E (E 11 , E 12 ) describe two more orthogonal modes (D y and Q x ), while the other elements (E 21 , E 22 ) specify two more orthogonal modes ( 9) www.nature.com/scientificreports/(D x and Q y ).When subjected to a 90-degree rotation, D x and Q y transform into D y and Q x , respectively.This establishes the polarization-independence of E modes.
To excite the dark resonances, it is necessary to reduce the symmetry of the UC.However, reducing the symmetry from C 4v to C 2v will not excite the dark modes.This conclusion is supported by Table S3 and Fig. S3 which display the IRREPs and electric field profiles of the C 2v group, and Table S4, which outlines the process of symmetry degeneration from C 4v to lower groups.However, it is possible to excite the dark resonances by reducing the symmetry of C 4v to C s .The IRREPs of C s are presented in Table 3.This reduction allows the existence of corresponding electric field components.In this case, according to Table S4, A 1 and B 1 (A 2 and B 2 ) of C 4v are reduced to A (B) of C s .
The symmetries of distinct non-degenerate modes correspond to different incident wave polarizations, resulting in polarization-dependent BIC modes even for the same perturbation type.As explained in Sect.S4 and shown in Fig. S4, BIC modes A 1 and B 1 of C 4v can be excited by an x-polarized plane wave, whereas BIC modes A 2 and B 2 can be excited by y-polarization, which indicates the possibility of selectively exciting BIC modes.Therefore, BICs are polarization-dependent.

Material property studies
Complex RIs of PEA 2 PbX 4 are shown in Fig. 2a-c.These materials exhibit neither dispersion nor loss in the THz region, making them suitable for employment in this spectral range.It is worth mentioning that the RI decreases as the halide changes from iodine to bromine and then to chlorine.This is because of the increase in the bandgaps (see Fig. S5).
Table 4 presents the mechanical properties of PEA 2 PbX 4 using the Voigt-Reuss-Hill approximations.The table shows that all three structures are mechanically stable (see Sect.S2) and ductile with B/G > 1.75 or ν > 0.26 according to DFT calculations.The bulk, shear, and Young's moduli exhibit an increase when the halogen atom changes from I to Br and then to Cl.This indicates that PEA 2 PbI 4 is more flexible compared to PEA 2 PbBr 4 and PEA 2 PbCl 4 .Notably, silicon and GaAs have very high Young's modulus values of 174.8 GPa and 87 GPa, respectively, making them stiffer than PEA 2 PbX 4 , which exhibits Young's modulus range of 12-22 GPa 41 .Therefore, PEA 2 PbX 4 -based metasurfaces can be used in wearable photonic devices.In addition, the relatively low bulk modulus of all three structures suggests their inherent softness, enabling easy transformation into thin films.This characteristic is particularly crucial for photonic applications.
Table 5 investigates the thermodynamic stability of PEA 2 PbX 4 by calculation of the formation enthalpy energy for its component materials.As can be seen, PEA 2 PbI 4 has a more negative value of formation enthalpy energy.This material is about 2 times more stable than PEA 2 PbBr 4 and PEA 2 PbCl 4 , respectively, and is some orders of magnitude more stable than 3D perovskites 42 .
Table 6 shows the adsorption energy of water molecules on PEA 2 PbX 4 (see Sect.S3).The proposed materials exhibit better hydrophobicity than most perovskites, making them more moisture-stable 43 .PEA 2 PbI 4 adsorbs water molecules approximately four times less than conventional 3D perovskites 44 .Among the proposed materials, PEA 2 PbI 4 exhibits the highest stability, whereas PEA 2 PbBr 4 shows the lowest stability against moisture.Therefore, adsorption energy and formation energy are directly correlated, and materials that have higher thermodynamic stability also tend to have better moisture stability 45 .
Based on the above results, we chose PEA 2 PbI 4 for the metastructure due to its higher RI in the THz region, superior thermodynamic stability compared to PEA 2 PbBr 4 and PEA 2 PbCl 4 , better moisture stability, and ductile nature.

Wave propagation studies and mode analysis
If d = 0, the structure exhibits C 4v symmetry, and its transmittance spectrum for x and y polarizations is depicted in Figs.3a and 4a, respectively.For both polarizations, two bright modes belonging to IRREP E are observed.Since q-BICs do not exist at d = 0, there is no energy leakage from the bound states to the zero-order channel.Furthermore, these modes are positioned at identical frequencies under both x and y polarizations, with their fields rotating 90 degrees from each other, highlighting their degeneracy (see Fig. S6).By adding an offset, it becomes possible to couple x-polarized waves with A 1 and B 1 , and y-polarized waves with A 2 and B 2 .This has been identified in the group theory analysis in "Structure design and symmetry analysis" section.As a consequence, the excitation of multiband q-BICs is selectively achieved depending on the polarization of the incident waves.In addition, symmetry breaking leads to the emergence of a transmission peak (mode A 1 ) within the near-zero transmittance valley of the non-BIC mode under x-polarization.The resonant frequency of A 1 is almost identical to that of the bright mode, producing an electromagnetic-induced transparency (EIT) effect.Moreover,  where ω 0 is the resonant frequency, a 1 , a 2 , and b are the constants, and γ is the total rate of damping that characterizes the Q-factor ( Q = ω 0 /2γ ) of the q-BICs.Figures 3b, c and 4b, c illustrate the fitting results for the four q-BICs at d = 0.5 µm.The Q-factors at d = 0.5 µm are 1.2 × 10 4 , 5 × 10 5 , 4.7 × 10 5 , and 4 × 10 4 for A 1 , B 1 , A 2 , and B 2 , respectively.The absence of clutter modes in the transmittance spectrum prevents interference with resonant states in the desired frequency range while maintaining high Q-factor q-BIC modes.Furthermore, all modes in this frequency range exhibit a spectral contrast ratio and modulation depth of approximately 100%.The spectral contrast ratio is defined as [(T on − T off )]/[(T on + T off )] × 100% and modulation depth is defined as [(T on − T off )/T on ] × 100%, where T off and T on represent the minimum and maximum transmittance, respectively 47 .These remarkable modulation depths and spectral contrast ratios significantly improve detection accuracy and switching efficiency in sensor and switch applications.www.nature.com/scientificreports/ To provide an intuitive representation of the resonances, we analyze the near-field distributions of the displacement current, electric, and magnetic field within the UC shown in Fig. 5.In mode A 1 , the head-to-tail configuration of magnetic moments is quite evident.The presence of current loops in the plane perpendicular to it provides clear indications of the magnetic toroidal dipole (MTD) mode (also see Fig. S7a).This mode can be considered as the combined result of both intra-UC and inter-UC moments.Non-parallel magnetic moments are evident in mode B 1 , which represents a magnetic quadrupole (MQ).In the case of A 2 , the field map reveals a clear vortex of displacement currents threading through the inner ring.Additionally, magnetic moments take the form of a vortex within the plane perpendicular to the current.This configuration enables easy recognition of an electric toroidal dipole (ETD) (also see Fig. S7b).Mode A 2 can also be considered as the combined result of both intra-UC and inter-UC moments.Finally, in B 2 , an electric quadrupole (EQ) mode can be identified by antiparallel electric moments.The fourth column in Fig. 5 displays a 3D representation of the electric field.The displacement current is represented by gray arrows, while the magnetic field is depicted by black arrows.
To perform a more comprehensive investigation of q-BIC properties, we utilize the Cartesian multipole decomposition technique (see Sect.S7) 48 .During this process, we integrate the displacement current density within the UC, to gain insight into the distribution of the electromagnetic source in the far-field.According to Fig. 6, the dominant multipole components in A 1 , B 1 , A 2 , and B 2 are MTD, MQ, ETD, and EQ, respectively.These results are completely consistent with Fig. 5.
The verification of the four q-BIC modes' symmetry-protected nature is supported by an eigenvalue analysis of a part of their band diagram.In this analysis, the k x component of the periodic boundary condition sweeps from 0 to π/p (Γ to X in the first Brillouin zone).As shown in Fig. 7a and b, while all four modes theoretically have ultra-high Q-factors at the Γ point in the first Brillouin zone, their values decrease sharply from the Γ point, demonstrating their symmetry protected-bound state in the continuum (sp-BIC) nature.Moreover, in Fig. 7a, we observe a high Q-factor at k x = 0.8π/p, which corresponds to off Γ-BIC and lies beyond the scope of this work.
The electric field profile of the modes at k x = 0 and k x = 0.5π/p can be observed in the insets of Fig. 7a and b.It is evident that the nature of each mode remains unchanged despite varying k x , however, the Q-factors differ.We now examine how the Q value of the q-BICs is influenced by the offset distance of the inner ring from the center, as illustrated in Fig. 8.The results reveal that the Q-factor of the q-BICs exhibit remarkably high sensitivity to changes in d.When d deviates from 0, we observe an exponential decrease in the Q-factor (Q ∝ d −2 ).For instance, at d = 0.5 µm, the Q-factor is significantly higher compared to that of d = 1 or − 1 µm, with a difference of nearly four orders of magnitude.This exponential decline in the Q-factor of the q-BICs agrees well with findings from previous studies where an increase in the asymmetry parameter resulted in a similar trend 49,50 .Consequently, manipulation of the Q value of q-BIC modes can be realized with asymmetric parameters.

Considerations of fabrication non-idealities
PEA 2 PbI 4 has extremely low absorption losses in the THz region.However, in the real fabrication process, surface roughness and defects can lead to absorption and scattering losses, which are major concerns when designing the metasurface.These losses are considered by adding the imaginary part of the PEA 2 PbI 4 refractive index (k) (i.e., the extinction coefficient) as n PEA2PbI4 = n-jk.Figure 9a and b and their insets show that for d = 0.5 µm as the k of PEA 2 PbI 4 rises, particularly when k approaches 10 -3 , there is a noticeable increase in the full width at half maximum (FWHM) of q-BICs and a reduction in transmittance intensity.However, the transmittance spectrum remains largely unaffected when k ≤ 10 −8 .Meanwhile, the FWHM of non-BICs does not change, and all resonance frequencies exhibit no shifts.The impact of meta-atom losses is more pronounced on reflectance curves compared to transmittance curves (see Fig. S8).The results in Fig. 9c and d show that substrate losses have a notably smaller impact than meta-atom losses.Consequently, the observed variations in Q-factor and transmittance intensity of q-BICs are relatively minor.
In the experiment, the actual Q-factor can be limited by the finite size of the array 51 .Periodicity among meta-atoms suppresses radiative loss through near-field coupling.To investigate the near-field coupling among the meta-atoms, we simulate how the resonance evolves for different array sizes (the number of UCs in the metastructure).The results for A 1 are presented in Fig. 10.The electric field magnitude rises in proportion to the size of the array.Moreover, the central rings of the array demonstrate a greater level of field compared to those located at the edges.When a 27 × 27 array size is employed, the Q-factor of metastructure is very close to that of an infinite array.To achieve a Q value comparable to that of an infinite structure, the minimum array sizes for modes B 1 , A 2 , and B 2 are 15 × 15, 17 × 17, and 25 × 25, respectively (not shown here).
To examine the impact of fabrication process tolerance on the q-BICs of the proposed structure, we conducted investigations into the Q-factor and resonance frequency of the modes within a ± 5% range of change in the structure parameters.Based on the findings presented in Table 7, it can be observed that when the parameters of the metastructure undergo a ± 5% change, the Q-factors of q-BICs remain in the same order.Consequently, in practical scenarios, one can expect to observe high Q-factors.To discuss the red and blue shift of the modes we can refer to the dielectric resonator's Mie resonance frequency explained as 52 : where V (x, y, z) is related to the size of the resonator, while μ and ε represent its permeability and permittivity, respectively.In addition, θ is a constant value that applies to a specific resonance.The rise (fall) in height, radius, and period can cause the device volume V (x, y, z) to increase (decrease), and as a result, the resonant frequency redshifts (blueshifts).Therefore, the shift in resonance frequencies caused by the aforementioned parameters is consistent with the Mie theory.

Conclusion
In summary, we use DFT calculations with the CASTEP module to obtain electrical, mechanical, and optical properties, as well as the moisture and thermodynamic stability of 2D perovskite PEA 2 PbX 4 (X = I, Br, and Cl).By comparing the results among PEA 2 PbI 4 , PEA 2 PbBr 4 , and PEA 2 PbCl 4 , we find that PEA 2 PbI 4 has a higher refractive index in the THz range.In addition, DFT calculations show that all three materials are ductile and flexible.The formation energy calculations indicate that PEA 2 PbI 4 is thermodynamically more stable due to its higher negative value.Moreover, the surface adsorption energy of water on perovskite reveals that PEA 2 PbI 4 exhibits higher hydrophobicity, making it more moisture stable than the other two materials.We propose and analyze the excitation of sp-BICs in a novel PEA 2 PbI 4 -based metasurface, where multiple Fano q-BICs with ultra-high Q-factors can be excited.We employ group theory to justify the excitation of BICs and their polarization dependence.Through eigenfrequency and frequency domain simulations, we validate the group theoretical analysis and investigate the near-field distribution as well as the far-field scattering of q-BICs.Based on the results, we can conclude that MTD-EIT-BIC and MQ-BIC can be excited by x-polarized, and ETD-BIC and EQ-BIC by y-polarized incident waves.The effects of meta-atom losses, substrate losses, the finite size of the array, and https://doi.org/10.1038/s41598-023-49224-9www.nature.com/scientificreports/ of q-BIC become broader as d moves away from zero.The typical Fano formula is utilized for fitting the Fano resonance curves of q-BICs 46 :

Figure 3 .
Figure 3. (a) The transmittance spectrum under x-polarization with various d, (b) Fano fitting of A 1 at d = 0.5, and (c) Fano fitting of B 1 at d = 0.5.

Figure 7 .
Figure 7.The Q values of (a) A 1 and B 1 , and (b) A 2 and B 2 along the path from Γ to X in the first Brillouin zone.

Figure 8 .Figure 9 .Figure 10 .
Figure 8. Q values of (a) A 1 , (b) B 1 , (c) A 2 , and (d) B 2 as a function of the inner ring's offset distance from the center.

Table 2 .
Eigenfrequencies of the proposed structure with C 4v symmetry.

Table 3 .
IRREPs of group C s .

Table 7 .
The effect of fabrication tolerance on resonance frequencies and Q-factor.